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Equilibration and GGE for hard wall boundary conditions 
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In this work we present an analysis of a quench for the repulsive Lieb-Liniger gas confined to 
a large box with hard wall boundary conditions. We study the time average of local correlation 
functions and show that both the quench action approach and the GGE formalism are applicable 
for the long time average of local correlation functions. We find that the time average of the system 
corresponds to an eigenstate of the Lieb-Liniger Hamiltonian and that this eigenstate is related to 
an eigenstate of a Lieb-Liniger Hamiltonian with periodic boundary conditions on an interval of 
twice the length and with twice as many particles (a doubled system). We further show that local 
operators with support far away from the boundaries of the hard wall have the same expectation 
values with respect to this eigenstate as corresponding operators for the doubled system. We present 
an example of a quench where the gas is initially confined in several moving traps and then released 
into a bigger container, an approximate description of the Newton cradle experiment. We calculate 
the time average of various correlation functions for long times after the quench. 


Introduction. Nonequilibrium many body physics is 
one of the most challenging areas of research of modern 
condensed matter physics. There have been spectacular 
advances in the field, driven by experimental studies of 
dynamics in optically trapped atomic gas systems, sys¬ 
tems with extremely weak coupling to the environment 
allowing a study of essentially Hamiltonian dynamics of 
time evolution m- Encouraged by these experimental 
advances there has been great theoretical activity in the 
area |10fED] , focused on questions like does a steady state 
emerge, how do local observables equilibrate, is there any 
principle which allows us to relate the steady state to the 
initial conditions? 

One of the most important recent experimental nuT] 
and theoretical results is that there is a relation 

between the initial state and the long time steady state 
for time evolution of integrable models. This was as¬ 
cribed to the fact that integrable models possess an in¬ 
finite family of local conserved charges {h}, in involu¬ 
tion, which include the Hamiltonian H, typically identi¬ 
fied with l 2 '. 

[H,h] = [h,h,]=0,H = l2, (1) 

These conserved quantities imply that there is a complete 
set of eigenstates for an integrable model which may be 
parametrized by sets of rapidities {ki} which are simul¬ 
taneous eigenstates of all A ■ For the Lieb-Liniger Hamil¬ 
tonian, the model which describes the Newton Gradle ex¬ 
periment [T] , the the action of the charges on these eigen¬ 
states given by: A |{fc}) = X) shown for 

the Lieb-Liniger gas, by following its actual time evolu¬ 
tion numerically and analytically, |22| , that at long times 
the gas reaches equilibration with the density matrix hav¬ 
ing no time dependence and becoming diagonal in the 
basis {ki}. 

How to describe this diagonal, time independent, den¬ 
sity matrix in general is an open question. It was pro¬ 
posed that the diagonal ensemble in turn |22l [24] takes 
the form of a generalized Gibbs ensemble (GGE) |24H36| . 

PGGE = ^ exp (- aj^ (2) 


with the ai, the generalized inverse temperatures, en¬ 
coding the initial state |4>o) through the requirement 

{Ii}final — {PGGeIi} — (^0 \^i\ ^o) — initial- ^ 1® 

a normalization constant insuring Tr [pgge] = 1- This 
interesting proposal, while valid for the case at hand, the 
repulsive Lieb-Liniger model, fails for models with bound 
states (string solutions) [37], a large class of models which 
encompasses, among others, the attractive Lieb-Liniger, 
the XXZ Heisenberg chain and the Hubbard model. 

When a GGE description is valid it provides an el¬ 
egant shortcut to the computation of correlation func¬ 
tions at long times, without having to explicitly follow 
the time evolution or to compute overlaps. Instead, hav¬ 
ing reached equilibration the correlation functions (of the 
Lieb-Liniger gas) at long times, or in this case the time 
average of the correlation functions, may be computed 
by taking their expectation value with respect to the 
GGE density matrix, e.g. (0 (t — )■ oo)) = Tr[pGGE'd]- 
It was further shown |3S] that the GGE ensemble is 
equivalent to an eigenstate, pgge — |{^o}) ({^o}|) for 
an appropriately chosen |{fco}) so that (0 (t —>■ cx))) = 
({fco}| 0 |{fco}). Another approach, the quench action 
approach |39| . is of more general validity but is more 
difficult to implement. It allows the computation of the 
diagonal density matrix in terms of the overlaps of eigen¬ 
states with the initial state, but such overlaps are hard 
to determine and are known only for few initial states. 
Again, it was shown that the resulting diagonal ensemble 
is equivalent to an eigenstate. 

Most of the work done on the Lieb-Liniger model was 
concerned with periodic boundary conditions, exceptions 
are |4D|. Real systems have finite extent with typ¬ 
ically a parabolic potential confining the particles. We 
will approximate this parabolic confining potential as a 
hard wall boundary. We will study the system in the 
limit where the system size L ^ oo, the number of parti¬ 
cles N scales with the system size, N/L = const, and for 
times much greater then the system size, t > L/vtyp {vtyp 
is a typical velocity). This regime is relevant for many 
experiments. We note that this time scale has been the¬ 
oretically considered for the Tonks gas before gniiiT]. 
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Time average. We shall consider circumstances where 
the system does not necessarily equilibrate in the long 
time limit and focus instead on the long time average of 
a local operator (observable) 0 evolving from the initial 
state I'h), 

( 0 )^ = i dt (^-l ^ ( 3 ) 

= ^ Ea I 0 1«) («I 

= Ea ('f'l A) (A|0|A) (A|vE-), 

we find it is given by a diagonal ensemble in the limit 
T —>■ 00 . Therefore the time averaged expectation values 
of local observables is given by the diagonal ensemble. 
Here |A) and |/c) are exact eigenstates of the system in 
the box. 

The system we shall study is the Lieb-Liniger Hamil¬ 
tonian describing the 1-D system of bosons with short 
range interactions [miiaiis]: 

L 

Hll = J dx^dxb^ {x)d^b{x) + c ( 6 '*' (j:) 6 (x))^| . (4) 
0 

Here ( x) is the bosonic creation operator at the point 
X and c is the coupling constant. Hard wall boundary 
conditions are imposed: 

V’ (a;i = 0 ,X 2 , ...xn) = 0 

i^{xi,X2,:XN = L) =0 (5) 

with ■)/) {xi, ...xn) the wave function of the bosons in the 
region a;i < 3:2 < ... < xn- 

The exact eigenstates of the Hamiltonian with the 
boundary conditions given in Eq. (§ are given by |42| : 

■0 (|fci|,.. IAatI) =^C'{£}0 (ei |fci|, ..£7 v |*w|), (6) 

{H 

where {e} corresponds to the 2^ sequences ej = ±1 and 

c (£1, ...en) = n n*<j (1 - e.ifcd+edfcd ) ’ 

^ (fci, .../cat) ='^A{P) xi < X2 < ... < Xn 

p 

with A (P) = ni<j (1 + kp^-kpj and the sum ex¬ 
tending over A^! permutations. These are the eigenstates 
with periodic boundary conditions. Furthermore the ra¬ 
pidities ki = ei\ki\ satisfy the Bethe ansatz equations 

m- 


kiL = 'Kni+^^ 
3¥=i 


arctan 



-I- arctan 


c 

ki -f kj 


These are exactly the same equations as for a doubled 
system of length 2L with twice as many particles hav¬ 
ing 2N with rapidities {si |fci|}. There is a one to one 
correspondence between states of a system with hard 


wall boundary conditions and states of a doubled sys¬ 
tem with periodic boundary conditions where all the ra¬ 
pidities come in pairs {k,—k} |35]. The Bethe Ansatz 
equations which determine the allowed rapidities {k} for 
the doubled system can be translated in a standard fash¬ 
ion |35] into a set of integral equations for the rapidities’ 
densities. We denote, for a given eigenstate |{A:}) of the 
doubled system, by Pp (k) the Bethe density of particles 
so that 2Lpp (k) dk is the number of particles in the in¬ 
terval [k, k + dk] of the doubled system. Similarly ph (k) 
denotes the hole density and pt (k) = pp (fc) -I- pk (fc) 
the total density. The number of states |{fc}), consis¬ 
tent with a given set of densities, {pp (k), ph (k)}, is 
measured by the Yang-Yang entropy |^, Syy ({p}) = 

+Pp(A:)ln(^)). The den¬ 
sities {pp {k) ,ph (k)} for the doubled system are deter¬ 
mined by the thermodynamic Bethe Ansatz equations 
which enforce the periodic boundary conditions: pt (k) = 
h + with K {k,q) = 

Time average, quench action and GGE action. The 
time average of a local observable, Eq.([^, can be rewrit¬ 
ten as jSnilll]: 

(0)t^oo = / ^ (^) {{k}\ 0 \{k}) 

(7) 

Here ({fc}|0|{A:}) is computed in the hard wall (non- 
doubled) system, and the quench action is given by: 

Squench {{p (fc)}) = J 9 "^ P) P (^) + {{p (fc)}) 

with Jg^{k)p{k) = ^log(|($ I {k})\), where g (k) = 
g{—k). The extra factor of 7 in front of S'yy ({p (P)}) 
comes from the fact that we are only considering states 
where the rapidities come in pairs {k, —k}. 

The time average of the Lieb-Liniger gas with hard 
wall boundary conditions corresponds therefore to a sin¬ 
gle eigenstate the one that maximizes the quench action 
|44| . Let us denote solution quasiparticle density for the 
doubled system as Pp (k) and the quasiparticle density for 
the original system as Pp (k) with Pp (fc) = 20 (fc) Pp (k) 
and p^ {k) = p^ (-fc). 

We proceed to convert the quench action into a GGE 
description of the system with hard wall boundary condi¬ 
tions and use it to compute time average of local observ¬ 
ables. We begin by determining its conserved charges. 
Since an eigenstate, see Eq. is given by a super¬ 
position of rapidities of the form {si j/cij} all the even 
conserved quantities satisfy the relation: 

/2n0 (Ifcil ,.. IfcATl) = 5] kfN i\ki \,.. |A:iv|) (8) 

Therefore {l 2 n} form a set of local integrals of motion and 
the quasiparticle density Pp (k) being symmetric in k is, 
in turn, uniquely determined by its even moments which 
correspond to its conserved quantities |45| in particular 
{l 2 n} are complete. Hence the even local integrals of 










motion, {l 2 n}, determine final state in terms of the GGE 
density operator, pgge = exp (-^ a 2 „/ 2 n)- The in¬ 
verse temperatures, a 2 ni found from the initial state |<i>) 
setting {l 2 n) = {l2n)GGE- local observable can be 

written in the form: 


(A) 



1 ^ 

1 ^ 



{^)gge = / ^ (^) ({fc}| 0 \{k}) (9) 

with Sgge = f p(A:) + ^InZ) -|- ^Syy- 

We can identify g (k) = Q; 2 nA^" + ^ (since both 

the quench action and the GGE are equivalent to a single 
eigenstate of the Lieb Liniger Hamiltonian (which corre¬ 
sponds to the extremum of the path integral in Eq. 
we establish that {Q)rp^ao ~ i^)GGE- conclude that 
the time average of Lieb-Liniger gas corresponds to a 
GGE density matrix where the conserved operators are 
the even local conserved densities. We further note that 
when considering an operator 0 with support far away 
from from the hard wall boundaries we may as well calcu¬ 
late {'d)GGE with respect to the doubled system. Indeed 
(®)ggb = Tr { 0 exp(-X;Q; 2 n.f 2 n)} for both systems. 
Since operators l 2 n are local, all correlation functions 0 
may be calculated by considering paths where the propa¬ 
gator is the quadratic piece of ^ cy. 2 nl 2 n while the inter¬ 
actions are given by the quartic and higher order pieces. 
We note that paths that cross the boundary of the sys¬ 
tem are exponentially suppressed when 0 is far from the 
boundary. 

Examples: 1. Newton’s cradle type - eigenstate ini¬ 
tial conditions. We will consider the following setup: 
there is a large trap of length L with hard wall bound¬ 
ary conditions in which there are multiple smaller traps 
of lengths Li moving with velocities Vi. Each of the 
smaller traps contains a Lieb-Liniger gas initialized in 
an eigenstate described by quasiparticle density {k) 
with Pp (k) = Pp {—k) (we note that thermal states also 
correspond to specific eigenstates |3S]). At time t = 0 
the smaller traps are turned off and the whole of the gas 
expands into the larger trap. We would like to find the 
quasiparticle density of the long time averaged final state. 
To do so we use the fact that all the even local conserved 
quantities are conserved during the quench, so we need 
to equate their values before and after the quench. We 
will show in the supplementary online information that 
in the thermodynamic limit we do not need to consider 
the edge effects for computing the local conserved quanti¬ 
ties. Therefore we need to find a symmetric quasiparticle 
density that satisfies the following set of equations: 



Figure 1: (A) The system is initialized in a state where 

two hard wall Lieb Liniger droplets of length Li moving with 
velocity ±V inside a large hard wall trap of length L. (B-G) 
The velocity distribution for the BEG bottom and the ground 
state quench top for a variety of quasiparticle densities and 
interaction strengths. The time average velocity distribution 
in red the initial velocity distribution before the quench is 
shown in blue. (B-D) V = 5, L = 1, kp = 1. (E-G) V = 5, 
L = 1 n = 1. The initial velocity distribution is computed 
in the supplementary online information and is shown in blue 
while the final velocity distribution is shown in red. the initial 
velocity distribution of the BEG is shown in the form of delta 
functions. Note collision narrowing in (B-D) and broadening 
in (E-G). 


bled system. A solution to this equation is given by: 



This solution allows for the calculation of various corre¬ 
lation functions for the system. Note that in the case of a 
periodic boundary condition we would have received the 
answer p^ (k) = ^ {k - ^Vi). 

Gonsider now the quench dynamics of a system consist¬ 
ing initially of two boxes of length I with N particles each 
in the ground state moving with of opposite velocities V 
and —V see Fig. HA). In experiment one typically mea¬ 
sures the probability distribution for the particle velocity. 
It is given by the Fourier transform of the field-field cor¬ 
relation function 

P{v,t) = j dxe~'‘^’’’ {x)b{0))^ 


L 


dkp^p (fc) 




ppt) 




The extra terms k }^Vi stem from the fact that under 
a boost to velocity V) the wave function is multiplied 
by eyip {I'Y^rriiViXi) with rui — We note that here 
Pp (fc) = Pp {—k) is the quasiparticle density of the dou- 


We will be interested in its time average. This exam¬ 
ple has many similarities to the experiment done by Ki- 
noshita et. al [T] where the system is placed in a parabolic 
confining potential and initialized in a state with some 
of the particles going left and some of the particles going 
right. Here we have replaced the parabolic confining po¬ 
tential with a hard wall box and do not therefore expect 
this probability distribution to match well with the one 
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measured in the experiment. The reason being that when 
confined by a harmonic potential the bosons move up and 
down the potential which slows them down and speeds 
them up periodically. In our setup the particles hit a hard 
wall and have their velocities reversed after the collision 
(as such they experience no intermediate velocities). As 
a result our calculation is expected to underestimate the 
probability of a particle having low velocity. 

We now proceed with the calculation: An impor¬ 
tant ingredient in calculating correlation functions is 
the occupation probability of the doubled box (fc) = 
. To calculate it we first calculate the quasi¬ 
particle distributions of the smaller boxes. The 
ground state total density pt of the smaller boxes, 
in the limit of large c is determined from |85| : 

27:Pt (k) = 1 + 1 f'Xp ( 9 ) leading to pt (fc) = 

9 {—kp, kp) ^ (1 + + o(^). Furthermore it is 

possible to obtain a relation between kp and N with 

kp = ^ +. which implies that the total 

particle density of the doubled box is given by: pp (k) = 
sir (1 + ^ (-^F + v,kp + v), therefore the 

final total density is: pt (k) = ^ (l + ^^5;) + •••■ oc¬ 
cupation probability: 


fL(k) = AL ^ e{-kp + v,kp + v) ( 10 ) 

v=±F 

with Ap = j (1 -I- (1 — We now proceed to 

compute the field-field correlation function (x )b(0)). 
We will only consider the case when the points x and 
0 are far away from the boundaries of the box so 
we may use the doubled system for all calculations. 
In terms of the occupation distribution, fpik), the 
correlation function is given by |46| : {x )6(0)) = 

/ if /l (k) (k) exp {-x f dufp (t) Pu {k)) + h.o.t. 

with w (fc) = exp (- ^ / dqK {k, q) f (q)) +. exp (- , 

and Fp = ^kpAp and K{k,q) = (p_gy 2_^_^2 = I- 

The generating function (fc) satisfies the equation: 
2fPu (fc) = - exp (- f fp (s) K {u, s) Ps (k)) - 1 

yielding for large c: (fc) = (l + exp (-^)) + 

exp(-?Ft) 

i — ^ (k-u) and Jfp{u)Pu{k) = 

_|T (1 + exp (-^)) + exp (-^). ^ Comb¬ 

ing all we obtain the velocity probably distribution: 


P 


[v) ^ ^ (- 1 )-^ arctanAjj(u) ( 11 ) 




with A±±(u) = Cp ((1 - Fp^S^=^^) (±^ ± kp) + 1 ) 
and Cp = 


27r 


AkfAL (^l+exp(^-^^;|T^^ 


Note that the velocity distribution Eq.(lll, see Fig. 
SB -D), underwent a collision narrowing. The distribu¬ 
tion is the leading order term for the set up of a hard wall 


trap. In a harmonic trap, as argued before, the probabil¬ 
ity of a particle having low velocity would be larger due 
to having to move up and down the harmonic confining 
potential. 

2. Newton’s cradle type, EEC initial conditions. A 
very similar scenario happens when we initialize the state 
in a collection of EEC’s each of length Li moving with 
velocity Vi inside a larger trap of length L. At t = 0 the 
smaller traps are released and interactions are turned on 
so that the system is described by a Lieb-Liniger Hamil¬ 
tonian with coupling constant c. The initial state is EEC 
and can be described by a quasiparticle density |44| : 


P 


BECi 

P 


{x) 


Ti^a{x,T) 

1 + a{x, t) 


( 12 ) 


Here a: = -, = — (where is the particle density) and 

a(x,T) = (4+^) Ji+2ix (4V^)- With J„ 

a modified Eessel function of the first kind of order n. 
Ey an argument similar to the one given above the final 
quasiparticle density is given by: 


pI (k) = 


E h 

2L 


ph^' 


k+-V, 


+ Ph^’ 


k--v. 


(13) 

More generally any translationally invariant quench 
that may be solved using periodic boundary conditions 
it is possible to define a box quench which may be solved 
analogously to Eq. (131 above. In the supplementary 
online information we show that for a quench with two 
boxes (each of length I with N particles in each box in a 
EEC state) with velocities V and —V inside of a box of 
total length L the velocity probably distribution is given 
by: 


P(u) ~ nBphlL^^x 

TT 

Hp Hp 

Hi + h^- VHp) ^ Hi+ \{v + VKp) 

Here, Kp = (l - , Gp = 2nBp 

+ + 

2n (1 - Gp+St^iW+l'j and n = f. See Fig. |^E-G). 

We note that the average velocity distribution has 
broadened as compared to its value at the start of the 
quench, while in the previous case ground state initial 
conditions the distribution underwent narrowing due to 
the collisions. 

Conclusions. We have studied a quench of the Lieb- 
Liniger gas on an interval with hard wall boundary con¬ 
ditions. We introduced a doubled system with periodic 
boundary conditions that is equivalent to the one on an 
interval. We have shown that the GGE formalism ap¬ 
plies to the computation of time averages of local oper¬ 
ators and that the even integrals of motion form a com¬ 
plete set of local conserved quantities. We have used this 
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approach to compute a quench where there are several 
small traps inside of a larger one and the smaller traps 
are released. We found that the quasiparticle density is 
additive. We have also calculated the expectation values 
of some local operators for this quench and in particular 


the time averaged velocity distribution. In the future it 
would be of interest to extend this work to models with 
bound states. 
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Supplementary online 
information 

Correlation functions (zero temperature case) 

We would like to calculate various cor¬ 
relation functions (6^ (0) (0) b (0) b (0)), 

(0) (0) (0) 6 (0) 6 (0) & (0)), and {p{x)p{0)) 

((6^ (a;) 6(0)) and the quasiparticle occupation probabil¬ 
ity /i (^) = %lk) already been partially done in the 


main text) for quenches discussed in the main text. We 
will consider the case where both x and 0 are far away 
from the boundaries of the box. In this case the problem 
becomes translationally invariant and all correlations 
may be calculated using a doubled box with doubled 
quasiparticle density (see the discussion below Eq. ([^ 
in the main text). As such we may use the results found 
in [151 - 150] . We will consider the initial conditions where 
the system starts with two small boxes of length I with 
N particles each, with each box cooled to the ground 
state. We will also assume that the boxes have velocities 
V and —V. To make the computations tractable we will 
work only in the limit of large c and only to leading 
order in -. 

C 

We now compute the local correlation functions using 
this occupation probability. We begin with the correla¬ 
tion function (6l (0) 6l (0) 6 (0) 6 (0)). It is given bv|47|: 


(6t (0) 6t (0) 6 (0) 6 (0)) = 21 ^ A (ki) J (^ 2 ) ^ 


{ki - k2f 


- ^ 2 ) + C 2 


s±Alkl('- 


(14) 


For convenience we will denote = z {^ + (l ~ z))- Where the last equality is in the limit kp ^ V c. 
Furthermore it is possible to obtain the density density density correlation function in the same limit. It is given by 

m- 


(6t (0)6t (0)6t (0)6 (0)6 (0)6(0)) ^6 J ^ J ^ J (fc,)(A)^)/,. (^ ^ 


(fci - k2f (ki - ksf (k2 - kzf 


^ ^ /|3 

8^3 ^ 


2kF\^ fv\^ 


(/ci - fe) -I-(fci - fca) +C^{k2-ks,) +c^ 

(15) 


Where again the last equality is true in the limit kp ^ 
F < c. 

We now repeat the calculation of the field-field corre¬ 
lation function (Jd^ (x) b (0)) (already partly given in the 
main text). It is given by |4t)| : 

/ elk 

^f^(k)e-^'^^co{k)x 
X exp (-xf dtfp (6) Pt {k) 


Here uj (k) = exp {-^ j dqK (k, q) f {q)) ^ 

exp(—^), where Fp = AkpAp and K{k,q) = 
-—= -. Furthermore the function P* (k) satisfied 
the equation: 


k — t + ic 


Using this expression it is possible to obtain that: 


P,(t)a-T(l + exp(-^ 


.exp (- 


2Fr 


nc 


{k-t) 


From this we obtain that / fp (t) Pt (fc) = 


2^Ptik)--j-^^exp^-Jfpis)K{t,s)Psik)j-l + ,A^exp(-?^). Comb 


(16) ing we obtain that : 
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Fl' 


{b^x)b{0)) 


27T 

2 exp(-f|; 


27r 


expl-^ (l + expf-^))) / /i(fc)e 


exp(-2Fj^/7rc) 


„ {-"t (l + “» {-"§)} I -^‘X 


X sin kpx { 1 — Fl 


exp (-2Fi/7rc) 


TTC 


I f, 77 exp(-2FL/7rc) 

cos —a; 1 - Fl - 

Z V TTC 


(17) 


We now proceed to the density density calculation. We know that the density density function is given by 

'ki — k2 + ic\ fp{ki,k2)' 


{P (a;) P ( 0 )) = P^ - ^ J dkifL (fci) w (fci) J dk 2 fL (fe) w (fca) y 

Here p (fci, fc 2 ) = —i (fci — k 2 ) + f dtfL (t) Ft (/ci, ^ 2 ). Here the function Ft (^ 1 ,^ 2 ) satisfies: 

ki — t + ic k 2 — t — ic 


exp {xp{ki,k 2 )) 


2TTFt (ki,k2) = 


-exp - / /i (s) F (s, t) Fs (fci, fc 2 ) - 1 


ki — t — ic k 2 — t + ic 
From this we obtain that Ft (fci, k 2 ) = — {ki — ^ 2 ) + .... Combing we obtain that 


(18) 


(p(a;)p(0)) =p^ - ^exp ^1 + J dfci/i (fci) J dk 2 fL (fc 2 ) ^1 - 2f ^ exp (^-ix (fci - ^ 2 ) ~ 


2 4 / 2Fl\ I . 2f, ft , Pl\\ 2fV f, . Fl\ 

— n^ -exp I-I -.‘ii’" ' ^ --111 


=P 


TTC J X^ 


1 f 2 Fl\ 
-exp I-1 




TTC J a;3 (1 + 


Fl\ 


TTC J 


sin kpx 1H-cos —x 1H- 


Fl\ 


TTC J 


1 ( 2Fl\ 1 . /, (V (.^FlW 

-(=xn I-1 — sm kpx 1 H-cos —x 1 H-x 


2 l / 2 

TT^^C V "KC J X‘^ 


TTC / 


TTC / 


x|ricos|x|l- 7 H (l+^H-cos(H| + lip) (l + g) 


+ 


2t, {cos (^x {1 - 7,) (^1 + ti)) +COS (^x {1 {1 + 


As such to leading order in 1/c we have calculated all the 
correlation functions for the two box quench. 


Correlation Functions (BEC) 

We would like to carry out similar calculations 
to the ones done above in the case when there are 
two boxes each of which is initialized in a BEC each 
of length I with N particles. The boxes are mov¬ 
ing with velocities V and —V (the container box is 
assumed to have size L). We will calculate the expec¬ 
tation values of the operators (61(0)61(0)6(0)6(0)), 
(61(0)61(0)61(0)6(0)6(0)6(0)), (6l(x)6(0)) and 

{p (x) p (0)). We will work in the limit of large c and to 
leading order in 1/c. We will also assume that both x 
and 0 are far away from the box boundaries so that we 
may use the doubled box system to do all calculations. 


The first step towards this calculation is to calculate the 
occupation probability of the BEC quench / (k) = ■ 

It is known that for large c the total quasiparticle density 
satisfies: 


Ilf 1 

Pt{k) = — + — Pp {q) dq = — + 
ZTT TTC J ZTT 


From this we obtain that 


f{k) = BLX p 


BEC 


k ^ 
^~2 


BEC 


2N 

ttcL 




(19) 


( 20 ) 


Here for future use we have defined Bl = y x 


L ‘ 

27r ^ ttcL 

Furthermore we note that for large c: p^^’^ {k) ~ 
F O (^) with n = ^ |l4]. Next we know that 

m-- 




















(21) 


(&t (0) &t (0) & (0) 5 (0)) = 21 ^ A (ki) J^fL {k2) ^ 


( fci - fca)^ 

- + c 2 


^ ^ 


Here we have assumed that n <C <C c. Furthermore we may calculate the density density density correlator 
similarly, it is given by m-- 


{b^ (0) (0) b^ (0) b (0) b (0) b (0)) ^ 6 


dki f dk2 f dk^ 


{ki - k2f 


(fci - ksf 


{k2 - ksY 


2tt / 27r 


^ fL{kl)fL{k2)fL{kY' ^^2,2(1, 1 ,^ 2 ,^,, ^ ^2 , 2 

[ki-k2) +c^(ki-k3) +c^{k2-k3) + 


^ -^B^rN 
87r3 ^ 




4 __c _^ 

4 n 2 2 n 

X - 


+ “ 1 ) ^ + 5?\/l + 5 ? + F “ (1 + ^ + ^) 

“ 1) ^ 5 ? (1 + :^ + n) 


(”) 


Here we have assumed that n ^ V c. We can now calculate the density density correlation function. We know 
that the density density function is given by |50| : 


(p(a;)p(0)) =p^ - 


1 

47r^ 


*i/l (fci) w (fci) / dfc 2 /L (fc 2 ) w (fc 2 ) 


fci - fc2 + ic\ f p (fci, fc2) 


ki — k 2 — ic J V ki — h 


exp(a;p(fci,fc 2 )) 


Here w (fc) = exp(—^ / dqK {k, q) f (q)) = exp(—with Gl = ^B^n. Here p(fci,fc 2 ) = —*(fci —fc 2 ) + 
/ (t) (fci) fc 2 )- Here the function Pt (fci, ^ 2 ) satisfies: 


27:Pt (fci,fc 2 ) = 


fci — t + ic k 2 — t — ic 
ki — t — ic k 2 — t + ic 


exp j Jl (s) if (s, t) Ps (fci, fc 2 )^ - 1 


( 22 ) 


From this we obtain that Pt(fci,fc 2 ) = (fci — fc 2 ) + Furthermore ~ —1 (l — 2^ (fci — fc 2 )) ~ 

— exp (—2); (fci — fc 2 ))- We now obtain that p (fci, fc 2 ) = —i (fci — fc 2 )(l+^). Combining we obtain that: 


(p(a;)p(0)) =p^ - 


2 + 2 cos(f(.(l+f^) + i)) 


47r^ 


7TC J \ 7TC 


Blx 


X J dkip^^^ (fci) J dk 2 p^^^ (fc 2 ) exp (fci - fc 2 ) + 


2 + 2 cos(^(x(l + ^) + g)) 


471^ 


exp-IH- 

V TTC J \ TTC 


P£ • exp —2n x 1 + 


Gl 


7TC 


(23) 


We would now like to calculate the field-field correlation function. It is given by 


(&'!' (x) b (0)) = J ^fL (fc) e (fc) X exp ^-x J dtfL {t) Pt (fc) 


Here w (fc) = exp J dqK {k,q) f (q)) ^ exp(-^), where and K {k,q) = = 

function Pt (fc) satisfied the equation: 

27:Pt (fc) = - ^_^^^% xp J Jl (s) K {t, s) Ps (fc)^ - 1 

Using this expression it is possible to obtain that: 


(24) 

Furthermore the 

(25) 


aWSi-^(l+exp(-?(^ 




TTC 


(fc-t) 
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From this we obtain that f /l (t) Pt (k) = —^ (l + exp (—^^)) + iGL^ exp Combing we obtain that : 


{b^ (a;) 6(0)) = 


exp 

f-G.) 

\ TTC / 


27r 

exp 

f-G.) 

\ TTC / 


TT 

X f 

exp 

f-G.) 

\ TTC / 


exp 


exp - 


Glx 

27r 

Glx 

27r 


1 + exp - 


1 + exp - 


2Gl\ 


■KC 

2Gl 

TTC 


J fL{k)e 

Bl cos 


— ikx ( 1 — Gl 


exp(-2Gj^ /ttc) 


l-G. 


exp {—2GlIt^c) 




exp 


Glx 


1 + exp - 


V 


X cos ( — ( 1 - Gl 


27r 

exp {—2GlI'^c) 


2Gl\ 


Blx 


nexp —2nx [1 — Gl 


exp {—2Gl/'^c) 


The velocity probably distribution is then: 


P (v)dxe *2^ (6^ (x) 6 (0)) ^ ni^L 


exp 


Hl 


Hl 


^HI + Kv-VKl) HI + Hv + VKl) 
Here Hl = ^ (l +exp (-^)) + 2n (l - and Kl=(i- 


(26) 


q-boson regularization 

We wish to show that the in the thermodynamic limit 
the edge contributions to the conserved quantities van¬ 
ish. To do so we need to introduce a q-boson regular¬ 
ization of the conserved charges EH EH- The q-boson 
system corresponds to M bosonic lattice sites with each 
site having operators i3„, Bj^ and Nn = that satisfy 
the relations B„Bl - q-^BlB„ = 1, 
and = S)). The q-boson Hamiltonian is given 

by: 

1 ^ 

= (27) 

n—1 

The system is integrable since the Hamiltonian may be 
derived from the following transfer matrix 

^=(c(A) D5)5)=r«(A).-i.(A) (28) 

With 

Here x = -\/l — q~^ and q = e^. There is an infinite 
family of conserved charges In, the first few densities cor¬ 
responding to these conserved charges are given by: 

(n) =x^BlBn+i 



(30) 


Furthermore the open q-boson chain is also integrable 
[53]. It is known that the Lieb-Liniger gas is a limiting 
case of the q-bosons, where the limit is taken as 6 —>■ 0, 
M5 = L, 7 = ^ and q = e^. We shall show that in the 
limit Li —> oo for any finite x Ibe edges give no contribu¬ 
tion to the conserved quantities. Indeed we notice that 
the conserved quantities are linear functions of the expec¬ 
tations of various operators e.g. i?„+i) with 

~ Li terms in the sum. Furthermore by translational in¬ 
variance each of the terms gives the same contribution 
e.g. 

r = ^ (j* (no)) — Boundary Terms (31) 

Here uq is some site in the middle of the q-boson chain. 
We notice that the expectation values of the boundary 
terms have absolutely no L dependence (they are just 
proportional to the expectation value of the density, den¬ 
sity density, field-field and related correlation functions 
which do not scale with L). Therefore in the limit that 
Li —)■ oo we have that 7® = ^ (J* (no)) and the bound¬ 
ary terms have disappeared. Since the Lieb-Liniger gas 
corresponds to a limit of the q-bosons we see that it the 
thermodynamic limit the boundary terms do not effect 
conserved quantities. 


Initial Correlations 

We would like to calculate the velocity probability dis- 

) tribution when the traps are initially released at time 
equal to zero. This would help us compare with the time 
averaged case. The experimentally accessible quantities 
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are most easily given in terms of an average velocity prob¬ 
ability: 


Pav{v) = ^ j dx J dye {b'^ (x) b (y)) (32) 

In the case of the BEC it is not too hard to see that 


lexp(-^) 
L 


Pav {v) ^ V^ (-!)■’arctanAij(?;) (34) 


27r 




Paviv) = ^n{5{v-V) + 5{v + V)) (33) 

In the case of the two boxes in their ground state, 
following a derivation given above we see that the velocity 
probably distribution: 


with A±±(u) = Cl ((1 - (±^ ± kp) + ' 

n. — 27r T. _ or,„ 


and Cl = 


results are used in Fig. Bb-g). 


, with Jl = 2kp. These 


to|e 





